%% Practical: Sensorspace OAT
%
% In this practical we will run through an analysis step by step, several
% of which require manual intervention. This will go through the following
% steps:
%
%   1) Prepare raw data for OAT analysis
%   2) Bandpass filter data and split into epochs
%   3) Compute a first level GLM analysis with OAT
%   4) Visualise results with FieldTrip
%
% We will work with a single subject's data from an
% emotional faces task (data courtesy of Susie Murphy) and perform an ERF
% analysis in sensor space. This dataset can be downloaded from:  
% 
% www.fmrib.ox.ac.uk/~woolrich/faces_subject1_data.tar.gz
% 
% Note that this contains the spm file:
% spm8_meg1.mat
% that is an SPM MEEG object that has continuous data that has already been
% SSS Maxfiltered and downsampled to 250 Hz. 
% 
% and
% espm8_meg1.mat
%
% which is an SPM MEEG object that has the same data epoched into the
% different task conditions.


%% SETUP THE MATLAB PATHS
%
% Sets the Matlab paths to include OSL. Change these paths so that they 
% correspond to the setup on your computer. You will also need to ensure 
% that fieldtrip and spm are not in your matlab path (as they are included
% within OSL).%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
% This cell sets the Matlab paths to include OSL. Change the osldir path so that it corresponds to the setup on your computer before running the cell. 
osldir = '/Users/andrew/Software/Matlab/osl/osl_full/osl1.5.0_beta';

addpath(osldir);
osl_startup(osldir);


%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
%
% This cell sets the directory that OAT will work in. Change the workingdir variable to correspond to the correct directory on your computer before running the cell.

workingdir = '/Users/andrew/Projects/OSL_test/preproc_manual/faces_subject1_data';

cmd = ['mkdir ' workingdir]; if ~exist(workingdir, 'dir'), unix(cmd); end % make dir to put the results in

clear spm_files_continuous spm_files_epoched;

%% SET UP THE LIST OF SUBJECTS FOR THE ANALYSIS
%
% Specify a list of the fif files, structual files (not applicable for this
% practical) and SPM files (which will be created). It is important to make
% sure that the order of these lists is consistent across sessions. 
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:

% fif_files{1}=[testdir '/fifs/sub1_face_sss.fif']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss.fif']; 
% etc...

% spm_files{1} = [workingdir '/sub1_face_sss.mat'];
% spm_files{2} = [workingdir '/sub2_face_sss.mat'];
% etc...
% set up a list of SPM MEEG object file names (we only have one here)
% set up a list of SPM MEEG object file names (we only have one here)
spm_files_continuous{1}=[workingdir '/spm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/espm8_meg1.mat'];


%% SETUP SENSOR SPACE OAT SOURCE RECON
% 
% This stage sets up the source reconstruction stage of an OAT analysis.
% The source_recon stage is always run even for a sensorspace analysis,
% though in these cases it simply prepares the data for subsequent
% analysis.
%
% In this example we define our input files (D_continuous and D_epoched) 
% and conditions before setting a time frequency window from -400ms before
% stimulus onset to +500ms and 1 to 100Hz. The source recon method is set to
% 'none' as we are performing a sensorspace analysis

oat=[];
oat.source_recon.D_continuous=spm_files_continuous;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[4 100]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.4];
oat.source_recon.method='none';

%% SETUP SENSOR SPACE OAT FIRST LEVEL GLM
%
% This cell defines the GLM parameters for the first level analysis.
% Critically this includes the design matrix (in Xsummary) and contrast
% matrix
%
% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:
%
% Each contrast is a vector containing a weight per condition defining how
% the condition parameter estimates are to be compared. Each vector will
% produce a different t-map across the sensors. Contrasts 1 and 2 describe
% positive correlations between the each sensors activity and the presence
% of a motorbike or face stimulus respectively. Contrast 3 tests whether
% each sensors activity is larger for faces than motorbikes.

Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';

oat.first_level.cope_type='cope';
oat.first_level.report.first_level_cons_to_do=[2 1 3];
oat.first_level.bc=[0 0 0];

%% SANITY-CHECK OAT SETTINGS
%
% As well as using the settings we specified in the previous cell, calling osl_check_oat  has filled in a bunch of other settings as well.

oat = osl_check_oat(oat);

oat
oat.source_recon
oat.first_level


%% RUN OAT
%
% The OAT structure you have created, oat, should be passed to the function osl_run_oat.m, to run the pipeline. 
% However, the OAT structure should contain a structure (oat.to_do), which is a list of binary values indicating which part of the pipeline is to be run. %
% E.g. oat.to_do=[1 1 0 0]; will run just the source recon and first level stages, whereas oat.to_do=[1 1 1 1]; will run all four (source-recon, first level, subject-level and group level). Here we are not doing any group analysis, so we need:  oat.to_do=[1 1 0 0];

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);


%% VIEW RESULTS
%
% It is highly recommended that you always inspect the oat.results.report, to ensure that OAT has run successfully. 
% Open the web page report indicated in oat.results.report in a web browser (there will also be a link to this available in the Matlab output). This displays diagnostic plots. 
% At the top of the file is a link to oat.results.logfile (a file containing the matlab output) - you should check this for any errors or unusual warnings.
% Then there will be a list of reports for each OAT stage. 
% Click on the "First level (epoched)" link to bring up the first level reports.
% This brings up a list of sessions. Here we have only preprocessed one session. Click on the "Session 1 report" link to bring up the diagnostic plots for session 1 and take a look.

disp('oat.results:');
disp(oat.results);
disp(oat.results.logfile);
disp(oat.results.report);

%%%%%%%%%%%%%%%%%%%%%%%%%
%% View the GLM design matrix 
% 
% Running the OAT analysis will create an OAT output directory (whose name is the name set in oat.source_recon.dirname with a “.oat” suffix added). This contains all you need to access the results of the analysis. It contains the analysis settings and the pointers to files containing the results for each of the pipeline stages that have been run. The oat can be loaded into Matlab with a call to osl_load_oat.
% oat=osl_load_oat(oat.source_recon.dirname);This will display an image of the GLM design matrix defined above.
% (NOTE that column 1 is motorbikes, columns 2-4 are faces)

% load first-level GLM result
stats1=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

figure;imagesc(stats1.x);title('GLM design matrix');xlabel('regressor no.');ylabel('trial no.');

%%%%%%%%%%%%%%%%%%%%%%%%%
%% VISUALISE USING FIELDTRIP
%
% Next we will use an osl wrapper around a Fieldtrip function to the
% results from contrast 3 (Faces>Motorbikes). This requires us to define
% several parameters in the S2 structure. Critically, the oat analysis and
% the number contrast within it.
%
% Note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGPLANAR'; % can also set this to 'MEGPLANAR'
S2.first_level_contrast=[3];

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, dats, fig_handle]=osl_stats_multiplotER(S2);
